Reduced dynamics of two oscillators collectively coupled to a thermal bath 
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We study the reduced dynamics of a pair of non-degenerate oscillators coupled collectively to a 
thermal bath. The model is related to the trilinear boson model where the idler mode is promoted 
to a field. Due to nonlinear coupling, the Markovian master equation for the pair of oscillators 
admits non-Gaussian equilibrium states, where the modes distribute according to the Bose-Einstein 
statistics. These states are metastable before the nonlinear coupling is taken over by linear coupling 
between the individual oscillators and the field. The Gibbs state for the individual modes lies in the 
subspace with infinite occupation quantum number. We present the time evolution of a few states 
to illustrate the behaviors of the system. 



I. INTRODUCTION 

We consider a composite system of two non-degenerate 
oscillators coupled collectively to a bath, i.e. the change 
in the occupation quantum number of one oscillator mode 
due to the environmental influence induces a correspond- 
ing change in the other mode. In the studies of environ- 
mental influence on a pair of oscillators, the oscillators 
are usually coupled separately to the bath through linear 
interactions [1]. In contrast, our main focus is on the 
three-body interactions between the oscillators and the 
field modes of the bath. 

The system is related to the trilinear boson model in 
quantum optical systems, where it is used to describe the 
process of parametric amplification and frequency con- 
version [2-5]. The two oscillators then play the roles of 
the pump mode and the signal mode, respectively, which 
are coupled to the idler or vibrational mode of a nonlinear 
medium in which nonlinear interactions are assumed to 
be dominant. By promoting the idler mode to a field and 
assuming that the coupling is weak, we can employ the 
standard open quantum system approach [6, 7] to study 
the damping of the system of oscillators in a thermal bath 
of the field. 

As far as we know the reduced dynamics of this model 
has not been discussed before for the general sectors in 
the non-integrable region. We find that the Markovian 
master equation as a time independent eigenvalue prob- 
lem can be solved analytically. Due to the collective cou- 
pling between the oscillators and the bath, the reduced 
dynamics exhibits the behaviors of finite-level systems 
[8-10] under the cascade process, even though we are 
dealing with a continuous variable system. 

The Hamiltonian of the system has the same for- 
mal structure as the Lee model for the bosonic systems 
[11, 12], and the one-particle sector of its oscillators sub- 
system is equivalent to the Friedrichs model [13]. The 
Friedrichs-Lee model was originally devised to study the 
effect of perturbation on the spectra in the Hilbert space 



[13], the mathematical structure of renormalizable quan- 
tum field theory [11, 14], and later on to study the non- 
integrable systems where resonance states emerge [15- 
17]. 

The reduced dynamics for this type of interaction 
shows a few unique features. It gives rise to a family of 
non-Gaussian equilibrium states confined to their respec- 
tive irreducible subspaces, whereas the Gibbs state of the 
individual oscillators is recovered in the subspace with 
unrestricted occupation quantum number. The oscilla- 
tor modes in the equilibrium eigenstates distribute ac- 
cording to the Bosc-Einstcin statistics [18]. These states 
are metastable before the nonlinear coupling is taken over 
by linear coupling between the individual oscillators and 
the field. 

In our discussion, we first present the Hamiltonian of 
the system in Sec. II. The Markovian master equation 
of the reduced system and its bosonic representation are 
then presented in Sec. III. We then solve for the equilib- 
rium states in Sec. IV, and study the time evolution of 
some states in Sec. V. 



II. THE HAMILTONIAN 

We consider a system of two oscillators and a field in 
one dimensional space, labeled by 1,2, and k, respec- 
tively. The free Hamiltonian is 



Ho = LUia\ai + uj 2 a\<i2 + 



^2w k a\a k ■ 
k 



(1) 
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where we use the units Ti = c = \, and uji,uj2 are the 
natural frequencies of the respective oscillators. We as- 
sume that wi > ll>2, and the field is massless, with the 
dispersion relation u>k — \k\. The creation and anni- 
hilation operators a\ and ctj obey the commutation re- 
lations, [ai,aj] = Si.j,i,j = 1,2, k. We normalize the 
field in a box with length f2, so that k = 2ni/Vl : with 
i = 1, 2, 3, • • • . The limit 51 — > oo will be taken eventu- 
ally, but we continue to use the discrete notation in the 
expressions below. 



We assume that the oscillators are coupled collectively 
to the field through the interaction 



v=xj: 






(L + a k + L- 



(2) 



in analogy to the linear coupling model between an os- 
cillator (labeled by a,al) and the field, i.e., a^au + aa\. 
A derivation of this interaction can be found in Ref. 4 in 
the context of trilinear boson model, where it is used to 
describe parametric amplification and frequency conver- 
sion in quantum optical systems [2-5]. The L± are the 
ladder operators of the SU(2) algebra, 



L A 



a[a2 



L_ = a\a 2 ■ (3) 

They raise and lower the r-quantum number of the com- 
posite system, respectively, see Eqs.(15)-(17b) below. In 
Eq.(2), A is a dimensionless coupling constant, and v(u>k) 
is a form factor that contains a high frequency cut-off to 
regularize the interactions. 

The total Hamiltonian of the system, H = Hq + V, 
when written explicitly in terms of the individual mode, 
has a form similar to the Lee model for bosonic system 
[11, 12]. Its one particle sector of the oscillators subsys- 
tem is equivalent to the Friedrichs model [13]. 

The system possesses two independent constants of 
motion, 



N = a[ai + a\a2 ■ 



N lk 



k 



(4) 



N remains a constant of motion of the reduced dynamics 
when the field modes are traced out. In the unstable 
regime, the system develops a resonance at the frequency 

ojq = oji - oj 2 , (5) 

where the 1-oscillator turns unstable and decays into the 
2-oscillator and the field [15]. Complex poles correspond- 
ing to the unstable oscillator then arise in the complex 
energy plane. 

From another point of view, the a^s can be regarded as 
the normal modes of two degenerate oscillators coupled 
through the SU(2) coupling interactions [19]. The sys- 
tem Hamiltonian H is then unitarily related to a family 
of Hamiltonians by SU(2) transformation. The details 
are presented in App. A. Therefore, the reduced dynam- 
ics discussed below is also applicable to these systems. 
It is then interesting to note that experimentally [20] it 
had been shown that spatial wave patterns generated by 
three-dimensional coherent waves obtained through the 
longitudinal and transverse coupling of laser modes in 
a cavity [21] are related to the eigenstates of the system 
with SU(2) coupling interactions. It may then be possible 
to realize the spatial profiles of a mixture of eigenstates 
of the system in the future. 

We could gain further insights into the conditions un- 
der which the model is applicable by comparing the 
Hamiltonian to anharmonic interactions [22]. However, 
since this comparison is outside the main line of our dis- 
cussion, we present it in App. B. 



III. MARKOVIAN MASTER EQUATION - 
BOSONIC REPRESENTATION OF SU(2) 

A. Markovian master equation 

The reduced dynamics of the two oscillators subsys- 
tem immersed in a thermal bath can be derived by trac- 
ing out the field degrees of freedom with standard meth- 
ods [6, 7], or through the complex spectral representation 
[23]. In the derivation, we assume that the oscillators and 
the fields are initially factorizable, and we use the weak 
coupling limit, or cquivalently, the A 2 £-approximation 
[22, 24] or the Born-Markov approximation [6, 7, 25]. 

Since the structure of the interaction Hamiltonian is 
similar to the linear coupling model between a single os- 
cillator and a field [7], the master equation acquires the 
Kossakowski-Lindblad (KL) form [26, 27], with a mod- 
ified unitary part. The reduced dynamics is therefore 
completely positive too [26, 27]. The reduced density op- 
erator of the two oscillators, /, evolves according to the 
equation df/dt = —Kf, where 



K = K + K d 



can be decomposed into a unitary part, 



(6) 



K f = i[H^,f], (7) 

Hq = (wi — 5bJi)a\ai + (0J2 — Suj2)a2 a 2 ~ 5bj' a\aia 2 cL2 
= (co - Slo )L - |K - 5u/ )N + Soj q (L 2 q - \N 2 ) , 

(8) 

and a dissipative part, 



K d f = -^n Q (2L + fL_ - L_L+f - fL_L+) 

- | 7 (f»o + 1)(2L_/L+ - L+Lj - fL + L_) . (9) 

In the first equality of Eq.(8), we have presented H' in 
terms of the creation and annihilation operators to ex- 
hibit the frequency renormalization SuJi to the individual 
oscillator. The explicit expressions of the coefficients are 



Sojq = Su>i + S0J2 . 
Sujq = Suji — 5uj2 ■ 
A 2 



8u>i — 



Z ^ / ■ 



5u 2 



Q/2ir * — ' uik — ^0 



l«K)| s 



i). 



Vl/2k 



Ep 



ujk — wq 



-Tlk- 



7 = 27rA 2 |«M| 2 



(10a) 
(10b) 

(10c) 

(lOd) 
(lOe) 



where ujq is the resonant frequency (5) and 7 is the decay 
constant. We note that the natural frequencies of the os- 
cillators are renormalized with opposite signs, compare 
Eqs.(lOc) with (lOd). In scattering problems, the aver- 
age number of field modes is zero n^ — 0. In this case, we 
find that only the frequency of the 1-oscillator is renor- 
malized, consistent with the discussion in Ref. 11. 



We assume that the field modes are in thermal equi- 
librium satisfying the Bosc-Einstein distribution n& = 
l/[cxp(ajfc/3) — 1], where (3 = l/(fcsT). For the resonant 
mode, we label its occupation number by 



1 



n 



e "o/3 _ i 



(11) 



B. SU(2) bosonic representation 

The generators of the SU(2) group in terms of the 
bosonic representation [28] are 



L Q = (a\ai - a\a 2 )/2. 
L 1 = (L+ + L_)/2, 



L 2 = (L + -L_)/2i. 



(12a) 

(12b) 



They obey the commutation relations [Li,Lj] — i^i^L^. 
The Casimir operator of the SU(2) is [29] 

L 2 = \{L + L_ + L_L + ) +L 2 = L+L_ + L (L - 1) , 

(13) 

which commutes with the LjS, [L 2 ,Li] = 0. The total 
occupation number of both oscillators as denoted by N 
(4) remains a constant of motion of the reduced dynam- 
ics. It commutes with the generator of the SU(2) group, 
[N, L^ = 0,i = 0, ±. The quantum number N will be 
used to label the irreducible representation, see Eq.(15) 
below. 

We make use of the occupation number basis 



„,„,.) = Mr MS„, 0) 



(14) 



and denote a state in the irreducible subspace labeled by 

N as 



r n 



\ni,n 2 ) 



where 



m - n 2 



N = m + n 2 . 



(15) 



(16) 



are related to the eigenvalues of L$ and L 2 through 
Eqs.(17c) and (17d) below, respectively. Using these la- 
bels we can establish the following relations, 



L+\r) N = \^{N + r + 2){N - r)\r + 2) N 
L-\r) N = y(N + r)(N-r + 2)\r - 2) N 

Lo\t)n = \t\t)n , 

L 2 \r) N = \N{N + 2)\r) 
N\r) N = N\r) N . 



IN 



(17a) 
(17b) 
(17c) 
(17d) 
(17e) 



Whenever a 1-oscillator is created, a 2-oscillator is annihi- 
lated, and vice versa. Consequently, the index r changes 
in step of ±2 under L±. There is a total number of 7V + 1 
substates in each irreducible subspace, and r ranges from 
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FIG. 1. Energy level diagram and the rate of transitions be- 
tween different levels in the p N > subspace (between diagonal 
elements), see Eq.(21) for the definition of fi ■ The transi- 
tion rates between two levels are proportional to the numerical 
constants on the left of the levels. Transitions between dif- 
ferent irreducible f^ N ' subspaces are forbidden by the SU(2) 
symmetry of the reduced dynamics. 



-N, -N+2, ..., N-2, N. The highest and lowest states 
are |iV)jv and \—N)n, annihilated by the raising and low- 
ering operators, L±\±N)n — 0, respectively. The state 
|r)jv has energy 



E = uirii + Lu 2 n 2 = \{Nlu' + ru>o) , 



(18) 



see Fig. 1 for a plot of the energy levels. They are non- 
degenerate if the ratio uji/lj 2 is not a rational number. 
We denote the basis in the Liouville space by 



f (N,N) 



l r )jvjv(^l = \ni,n 2 )(mi,m 2 \ 



(19) 



We find that using this notation is more convenient for 
our later discussion since it is more compact and it man- 
ifests the fact that N and N are constants of motion 
of the reduced dynamics. Since the L± operators come 

in pairs in Kd (9), Kf^.~' is a linear combination of 
f;.f ,fr±2-?±2- Consequently, the quantity 



(20) 



is a constant of motion under K, and the basis states in 
each f( N > N > subspace arc connected to the others with 
the same v value only, under the reduced dynamics. 

It is interesting to note that the generator of the time 
evolution K is invariant under a rotation along the L$ 
axis, as shown in App. D. Furthermore, the SU(2) gen- 
eralized coherent states [30-32] reside in the f( N > N ) sub- 
space, see App. E for details. 



IV. EQUILIBRIUM STATES 

In this section, we obtain the set of equilibrium states 
for the irreducible subspaces. These states are non- 
Gaussian states in the coordinate space, as opposed to 
the Gaussian Gibbs state for an oscillator in a thermal 
bath. We will also show that the oscillator modes in 
these equilibrium states distribute according to the Bosc- 
Einstcin statistics, and the Gibbs states of the individual 
oscillator are recovered when the occupation quantum 
numbers for these modes are not restricted. 

We first note that the equilibrium states have non-zero 
trace and can be decomposed into the diagonal basis state 
labeled by 



/, 



(N) = 



■(N,N) 
r;r 



(21) 



whenever N — N and f — r. We can then write the 
equilibrium state as a sum of the diagonal basis states, 



AN) 

</cq 



l N 

— ^PN-Inf, 

■>N i — * 



(N) 

N-2n ' 



(22) 



where Zn is a normalization constant, and pn-2u are 
coefficients to be obtained below. Note that since each 
fr = |ni)(ni|®|n2)(n2| is separable, /i q is a separable 
state. 

The action of K on these states is 



KfW = u r f%l + v r fW+w r ri N _ 



(N) 
2 ) 



where 



= -±no[(N + l) 2 -(r + lf], 



(23) 



(24a) 



l(2n + l)[(N + i) 2 -r 2 -l]+lr, (24b) 

(24c) 



-l(n + l)[(N+lf-(r-lf 



by using Eqs.(17a)-(17e). The transition rates between 
different energy levels within the same f^ N ' subspace 
are summarized in Fig. 1. They can be read off from 

Eqs.(23)-(24c) directly to give 



two-level systems collectively coupled to a field [33, 34], 
supcrradiance [33] may also arise in this model. How- 
ever, the energy levels for the two models are different, 
compare Fig. 1 in this paper to Fig. 1 in Ref. 33. This 
is because a symmetrization of the underlying states is 
required for N two-level systems but not for a pair of 
non-degenerate oscillators. 



The equilibrium condition, Kf cq 
as a matrix equation, 



(N) 



K-fW=0, 



0, can be written 



(26) 



where 



K 



/V N U N - 2 

W N VN~2 Wat_ 4 

W N ^ 2 «JV-4 WiV-6 



W 6 -N V4-N U2-N 







W 4 -N V 2 -N U-N 

w 2 -n V-n/ 

(27) 



is a tridiagonal square matrix of dimension N + 1, and 
the eigenvector is a column matrix, 



/ Pn \ 

PN-2 
PN-4 






1 
Zn 



(28) 



P4-N 

P2-N 

V P-N J 



The solution to the coefficients is 
PN-2n 2 -n^-" 2 (l + n )" 2 , 712 = 0,1,2, 



,N, 
(29) 



which can be checked to satisfy Eq.(26) easily. By requir- 
ing Tr(/ e q ) = 1, we obtain the normalization constant 



/, 



(TV) 



(JV) _ I , if r = N. 



Jr + 2 \ u r , if r^N 

-> / W = Vr , all r 

AN) _ j w r , if r ^ -N , 

Jr ~ 2 if r = -N. 



(25a) 

(25b) 
(25c) 



The maximum rate is proportional to the numerical fac- 
tor N(N + 2)/4 for a transition between the /± 2 and 

the /q levels for even N, whereas the maximum rate is 
proportional to (N + l) 2 /4 for a transition between the 

/ x (Ar) and the /[^ levels for odd N [7]. 

Since the transitions between the energy levels in the 
f( N ^ subspace are very similar to the Dicke model for N 



N 

?N = / J PN-2n 2 — ! ' "T "() i - "(| 

n 2 =0,l,2,--- 



(l + rl ) jV + 1 -n Ar+1 , (30) 



r N ) = 



by using the relation (1 — r)(l + r + r + ■ ■ 

I _ r N+l 

Furthermore, the equilibrium occupation probability 
P2n 1 -N can be cast into the form (for fixed N), 



,N 



PAT-2„ 2 = [e^ 2 (l+n )] e- fjE (31a) 

(31b) 



OC g _ ™l' 3i * J l e _n 2/3w2 



where E is the energy of the configuration (18). Eq.(31a) 
shows that pN-2n 2 satisfies the canonical distribution for 



fixed N, whereas Eq.(31b) shows that the system obeys 
the Bosc-Einstein statistics with two energy modes [18], 
where the 1- and 2-oscillator states play the roles similar 
to the excited and ground states in a two-level system, 
respectively. 

A special case occurs at zero T, or equivalently, when 
the field mode has zero occupation number, no = 0. 
Damping to the system can then be caused only by the 
spontaneous emission of a field quantum [3] , accompanied 
by the lowering and raising of the 1- and 2-oscillator's 
occupation quantum numbers, respectively. In this situ- 
ation, the lowest state in each subspace, such as f_ N in 

the diagonal subspaces, and p : '- for N ^ N in the 

off-diagonal subspaces, are annihilated by Kd- Hence, 
the subspace spanned by these states does not experi- 
ence decoherence. Therefore, they should be included in 
the expression of the equilibrium state for T — 0, 



•'cq 2-^i 






AT=0,1,2, 



NjtN 

E 

AT,JV=0.1,2. 



(N,N) 
NN- 1 -N.-N 



if 



(32) 



where cn are real coefficients subjected to the normal- 
ization condition ^2 N c/v = 1, and c' ~ are complex co- 
efficients constrained by the positivity condition of /' . 
/' is a separable state by inspection. Indeed, since 



/ 



(N,N) 
-N,-N 



|0,JV)(0,JV| = |0)(0|®|7V)(JV| 



(33) 



we can factor out the overall 1-oscillator state |0)(0| from 
the right hand side of Eq.(32), so that /^ q = |0)(0| 
p( 2 \ where p^ is the density matrix of the 2-oscillator. 
Hence, it is clear that f is a separable state. 

It is interesting to compare the reduced dynamics of 
the two oscillators system to that of a single oscillator 
coupled to a bath. In the case of a single oscillator, the re- 
duced dynamics too has the Kossakowski-Lindblad form 
[26, 27], although it has a different set of L\ operators, 
L = a^a + 1/2, L' + = a', and L'_ = a, where a\a are 
the creation and annihilation operators of the oscillator. 
When written in terms of super-operators, the reduced 
dynamics has the SU(1,1) symmetry [35]. This reduced 
dynamics connects all the number bases in this system, 
and the equilibrium state is the Gibbs state. 

For the system we consider, if we trace out the 2- 
oscillator state from Eq.(22), we find that the 1-oscillator 
equilibrium state becomes 



AN) = (l + no) N 

•/cq.l 



'N 



N 

E 

m =0,1,2, 



e- ni *""|ni;ni». (34) 



When N is unrestricted, we recover the Gibbs state. 
Tracing out the 1-oscillator state will produce a simi- 
lar expression of / eq I as in Eq.(34), except \ni\n\f) is 
replaced by \N — n<i\ N — n 2 )), whereas the rest of the 
coefficients remain unchanged. 



V. TIME EVOLUTION OF STATES 

To study the time evolution of the reduced system, 
we first introduce the interaction picture for the reduced 
dynamics. We denote the density state in this picture as 



p = e Kat p = e^Ve- iH o* . 



(35) 



By expanding p = J^ c,/j and p — ^2 i c,/, in terms of 

the time- independent basis fi = /,:.-' , the coefficient c, 
acquires a phase due to the action of exp(iKot) on /j, 

5j = c i exp(z6'it) , (36) 

in which 

Oi = \H- 5co' )(N -N)+ i(« - 5iu Q )(r - r) 

-\5u' {N 2 -N 2 ) + \8u:' {r 2 -? 2 ). (37) 

Since the phase angle vanishes for the coefficients asso- 
ciated to the probability elements fr , we have 5j = q. 
In this situation, we drop the tilde sign on the coefficients 
to simplify the notation. 

In the interaction picture, the equation of motion be- 
comes 



d 



-KdP, 



(38) 



where the effect of Kd on the basis state is 

iScu'r.i't ; 



K d fi 



■hn {2e l5 ^ vt L+f i L_ - L_L + f t - / 4 L_L H 



- i 7 («o + l){2e- i5 ^ vt L_f i L + - L+L_fi - ftL + L_) , 

(39) 

in which v is already defined in Eq.(20), see App. F for 
details. The extra phase factors in front of two of the 
terms L + /j£_ and L_/jL + are due to the action of L\ 
in Hq of Eq.(8). For basis states that lie in the probability 
subspace, we have v = 0. In this case, Kdfi reduces to 

K dfi- 

By noting that the /( JV > JV ) and f( N < N ) subspaces are 



related by the following relations, f~. r 

Kofff' N) = [K f$">V, and K d f^ 
we obtain 



(N,N) 






(iV,JV) n 



c(N,N) 



Kf^> = [Kf^>] 



(AT.ATht 



(40) 



Hence, we can deduce the action of K on one subspace 
from the other. We will next illustrate the general fea- 
tures of the time evolution of the system with a few ex- 
amples. 



A. States in lower subspaces 

We begin by studying the time evolution of a state 
dwelling in the subspaces up to N, N = 1. The density 



matrix is 



Pit) = Wo 



(0) 



•<*(*)/} 



(i) 






mm 



h.c. 



mn 



(1,0) 



h.c. 



(41) 



where the coefficients are subjected to the normalization 
condition d{t) + a(t) + b(t) = 1 and the positivity con- 
ditions of p. The effect K^fi can be worked out using 
Eqs.(39), (17a) and (17b). By equating the coefficients 
associated to the same basis state /, on both sides of 
Eq.(38), we find that the coefficients evolve as 



d = 0, 

a = -7(1 + n )a + -fn a b . 
b = 7(1 + n )a - jn b , 
c = -7(n + \)c. 
g= -h(n + l)g, 



h 



hjnoh, 



(42a) 

(42b) 

(42c) 
(42d) 
(42e) 

(42f) 



where we have omitted the time dependence on the coef- 
ficients for simplicity. 

We first make a few observations. We find that only 
the coefficients under the same f( N > N ) subspace are con- 
nected. We also find that the coefficients for the ft 1 ' 1 *) 
subspace, namely, a, b, c, evolve in exactly the same way 
as the amplitude damping channel for qubits [8] . In gen- 
eral, it can be shown that the f( N > N > subspace evolves 
similar to the TV-level system under the cascade process 
with a single decay constant in vacuum, 7, 



I TV) «• \N -1) 



o 



**\1)+* 



(43) 



as depicted in Fig. 1. For instance, the dynamics in the 
j( 2 ,2) su b S p ace behaves similar to the three-level system 
[10] under the cascade process. 
The solutions to Eqs.(42a)-(42f) are 



d(t): 


= do, 







14a) 


a(t)-. 


1 + 2n 


1- 


_ e - 7 (l+2fio)4" 


14b) 


b(t): 


= 1 - d Q - a(t) , 


(44c) 


m- 


= g oe -(2fto+i)7*/2 f 


(44d) 


~9(t)- 


= 5oe -7(«.+l)t/2 ) 


(44e) 


h(t): 


= h a e'^ ot/2 , 




( 


44f) 



where do denotes the value of d at t — 0, and etc. From 
these expressions, we learn that the diagonal compo- 
nents eventually settle down at some equilibrium values, 
whereas all the off-diagonal components vanish asymp- 
totically. In the general /( JV > JV ) subspace, we find that the 
diagonal coefficients contain the time exponential factor 
cxp[— Nj(2fit) + l)t], hence they evolve towards the equi- 
librium value more rapidly, whereas all the off-diagonal 



coefficients vanish asymptotically. However, as already 
shown in Sec. IV, for the special case of zero tempera- 
ture, the lowest off-diagonal coefficient in each subspace 
does not experience decoherence, as is clear from the ex- 
pression of h(t) in Eq.(44f), which is independent of time 
when no = 0. 



B. States involving infinite number of subspaces 

We now investigate the time evolution of states involv- 
ing the general subspaces. In principle, the initial states 
can be decomposed into basis states in the various f( N > N > 
subspaces, and the time evolution can be analyzed sub- 
sequently. Since this is a tedious and not an illuminating 
process, we will discuss the general features of the time 
evolution by comparing the initial and the equilibrium 
states with two examples. To simplify the expressions, 
we will make use of the dimensionless position coordi- 
nate 



rriiUJi 



(45) 



where rrii is the mass of the i-oscillator and ^ is the or- 
dinary position coordinate with the dimension of length. 
(1) Assume that the 1-oscillator is initially a super- 
position of two Gaussian states centered at x± = ±a, 
respectively [36], 



ii(a;i) = iV ie -( Xl -°) 2 + 7V 2 e-( Xl + a ) 2 = V c n ( Xl \n} 



n=0 



(46) 



where N\ , -/V2 are the normalization constants that 
give the relative height between the two Gaus- 
sians. In Eq.(46), we decompose 4>\ in terms 
of the harmonic oscillator wave function (x\ri) = 



H n {x) exp(— x 2 /2)/y / 2 n n!v / 7r, with the expansion coeffi- 
cient c„ = (n\4>i), and H n (x) is the Hermite polynomial. 
As an example, we choose a — 2 and N1/N2 — 2. A little 
calculation shows that ci is the dominant term, 

{c , ci, c 2 , • • • } = {0.34, 0.22, 0.78, 0.23, 0.41, 0.05, • • • } . 

(47) 

The density matrix p\ [ ni (xi,Xi) = (xi\4>i) (4>i\xi) in the 
coordinate space is plotted in Fig. 2(a). 

Consider an initially uncorrelated composite system of 
two oscillators, $ ini , with the 2-oscillator initially in the 
ground state, p^ ni = |0)(0|, 



$ 



(i) -_ (i) 

- Pl,ir 



P2.ini 



n, m— 



/ J c nCrnJi 



(n.m) 



(48) 



n,m— 



(e)p^ 




mation carried by the 1-oscillator in the c^s is inherited 
by the 2-oscillator to some extent, as can be seen in the 
density matrix of the 2-oscillator, 



P (1) 

t J 2,cq 



TrM 1 ! =J2\cn\ 2 W)(N\ 



(50) 



N=0 



in which only the highest energy level in each of the 2- 
oscillator subspaces is occupied. Since ci is the dominant 
term, p 2 w |c2| 2 |2)(2|. The plot in the coordinate space 
then gives a characteristic three-peak profile of the wave 
function (x\2) along the diagonal, see Fig. 2(f). 

When temperature increases, the populations of the 
1- and 2-oscillator start to distribute accordingly among 
different levels in the f( N ~) subspace. The reduced states 
are then approximately given by 



o (1) 

rl,eq 



o (1) 

P2,eq 



l£2| 

z 2 

N 
z. 



.( P ( 2) |2)(2|+ P ( 2) |l)(l|+^|0)(0|), (51) 
.(p^|2)(2|+^ 2) |l)(l|+ p ( 2) |0)(0|). (52) 



As a result, the three-peak profile is smoothed out, as 
shown in Figs. 2(c) and 2(g) for no = 1. By comparing 
Eq.(51) with (52), we notice that the order of the level 

(2) (2) 

populations between the 1- and 2-oscillator, i.e., p 2 ,Pq 

(2) 

and p_ 2 , i s reversed. This feature recurs in the other 



subspaces as well. For large temperature, p. 

for all i. Consequently, we obtain a uniform distribution 



(JV) 
i,eq 



5 N 



(1) 



eq 



(h 



(i) 



cq- 



as can 



among all the levels. In this limit, p\ 
be seen by comparing Figs. 2(d) and 2(h) for no = 10. 
(2) We next consider an initially entangled state, 



FIG. 2. Initial and equilibrium configurations of the indi- 
vidual modes of example (1) for several temperatures, with 
parameters a — 2 and TV1/TV2 = 2 in Eq.(46). 



The spatial profile of p 2 ini (xi,X2) is plotted in Fig. 2(c). 
Under the reduced dynamics, the off-diagonal coefficients 
undergo exponential decay, and $; ni eventually evolves 
into the equilibrium state 



*$ - E l<*l a /£ 



(JV) 



JV=0 



JV 



E m 2 E 



PN-2n 2 



7V=0 



n 2 =0 



^A' 



(49a) 



| TV - n 2 ,n 2 )(N - n 2l n 2 \ 



(49b) 



cf. Eq.(22) for the expression of / Cl 



/ oq 



In the special case of zero temperature we have no = 0. 
The 1-oscillator then settles down to the ground state 



p\ ' = Tr 2 $eq = |0)(0| with a Gaussian profile as de- 
picted in Fig. 2(b). Notice that the off-diagonal peaks 
in Fig. 2(a) have decohered away, much like its single 
oscillator counterpart in a thermal bath [36]. The infor- 



$. (2) ) = 



n=0 



Cn\n,n) 



(53) 



For a comparison with the first example, we choose c„ to 
be the same as those in Eq.(47). In this example, both 
the oscillators are initially in the same state, 



(2) 
/Vim 



(2) 



E 

n=0 



|c„r|n>(n|«M 2 |2>(2|. (54) 



The density matrix of the system, 
evolves into the equilibrium state 



^ ] = |$^)($. (2) | 

mi I mi / \ lni I 



*S> = E i c ^i 2 /< 



(2JV) 
cq 



(55) 



N=0 



which has been shown to be separable in Sec. IV. Hence, 
the entanglement between the pair of oscillators are lost 
eventually. In fact, any initial entanglement in the re- 
duced system vanishes asymptotically in view of the sep- 
arability of the equilibrium state / cq (22). 

A comparison of Eq.(55) with Eq.(49a) also shows that 
the choice of the superpositions \n,n) in Eq.(53) have re- 
sulted in the exclusion of the odd f( 2N+1 ) subspace from 



the equilibrium state. This is a consequence of the fact 
that the set of odd number subspaces are absent from 
the initial state. Therefore, by specifically preparing the 
initial state, some subspaces could be excluded from the 
equilibrium state. In this model, states with different ini- 
tial conditions may evolve into different classes of equi- 
librium states. 



VI. CONCLUSION 

We have considered a system of two oscillators collec- 
tively coupled to a field through a three-body interac- 
tion. The model is applicable within a time frame in 
which nonlinear interaction is dominant. The two os- 
cillator modes become effectively coupled as a result of 
their collective interaction with the field. Consequently, 
the reduced dynamics possesses the SU(2) symmetry 
that is common to finite-level systems, which leads to 
non-Gaussian equilibrium states for the collective modes. 
These are metastable states until linear interactions be- 
tween the individual oscillators and the field takes over 
and drives them to new equilibrium states. The re- 
sults suggest that new forms of equilibrium states could 
emerge when subsystems are collectively coupled to the 
environment under different symmetry of the reduced dy- 
namcis. 

It is interesting to further explore the implications and 
manifestations of the results in other systems, such as 
in the orbital motion of two-electron quantum dots [37] 
and light-phonon systems [38] . In the later case, the ex- 
istence of metastable states of photons may prevent pho- 
tons from thermalizing too rapidly with phonon bath. 
This may have interesting implications in photosynthetic 
systems, where the role of phonon is played by the vibra- 
tional modes of photosynthetic reaction centers [39] . We 
leave these interesting investigations to future works. 
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Appendix A: Oscillators coupled by SU(2) coupling 
interactions 



Using the operator [19, 20] 

U = exp(— i4>L ) exp(— i9L 2 ) 



(Al) 



where Li are the generators of the SU(2) group (12a)- 
(12b), the ais are related to a corresponding set of bi 
operators by 



«i 

02 



U 



r/t 



e 4 ^ 2 cos(0/2) e-^/ 2 sin(0/2) 
-e^/ 2 sin(0/2) e^ 1 ^ 2 cos(0/2) 



(A2) 



The family of Hamiltonian related to the oscillators sub- 
system is then given by 



(A3) 



#12 = wiai T ai + W2«2«2 

= Y(&i t &i + fr 2 t &2)+w o cos0£o 
+ ujo sin 0(cos <pL\ + sin <$>Lq) , 

with the corresponding interaction 

V ~ a\a 2 ak + a\a\a\ 

= [ — sin0L o + cos0(cos(/)L 1 + sin(j)L2)] (a^. + a* k ) 
- i (sin (pLi - cos (j>L 2 ) (au - a\) , (A4) 

where 0.1,0,2 in Li (12a)- (12b) are replaced by 61,62, re- 
spectively. The second equality is obtained by substitut- 
ing Eq.(A2) into the first line of Eq.(A4). Note that TV 
(4) also commutes with U. Hence, the total occupation 
quantum number of the oscillators remain the same un- 
der the change of basis. Note also that the vacuum state 
does not change under U, i.e., |0,0)' = C/t|0,0) = |0,0). 
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Appendix B: Comparison with anharmonic 
interactions 

The anharmonic interactions couple the position oper- 
ators of the three species. Using dimcnsionless position 
coordinate defined in Eq.(45), the position operator Xi is 
related to the a^ by Xi — (a^ + aJ)/V2- The anharmonic 
interactions then take the form 

V = \y^=$=( ai + a\)(a 2 + <4)(o fc + 4) . (Bl) 



k 



VW^f 



In this situation, we find that aside from the uoq resonant 
mode, there is a second resonant mode that occurs at a 
higher frequency to' compared to ujq, 



Wn 



Wl +L0 2 



(B2) 



This mode is initiated by the interaction terms a\a2<i\ + 
h.c, where h.c. denotes hermitian conjugate. If the occu- 
pation number of the uj' mode in the field is small enough 
then this mode is not excited. This could be achieved if 
the temperature of the bath, T, is low enough so that the 
condition 



UJ 



< k B T < w' Q 



(B3) 



is satisfied, where ks is the Boltzmann constant. An- 
other way to achieve this is to formally impose a high 
frequency cut-off w c to the form factor v(cok) in the in- 
teraction (Bl) so that 



k B T < w c < Wq 



(B4) 



is satisfied. Both conditions are separately consistent 
with the A 2 £-approximation [22, 24], or the Born-Markov 
approximation [7, 25] , used to derive the Markovian mas- 
ter equation of the system. Condition (B3) implies that 
Wi, W2 should be of the same order, since u>i — LU2 < k B T 
for small T, whereas condition (B4) implies that they 
should not be too small, since we require uj c <C Wi + ui2 
for a large cut-off frequency. In the latter situation, since 
wi and u>2 are not small, only a few lower energy modes 
will be excited in the dynamics, see Fig. 1 for the energy 
levels of the system. 

Aside from the resonant modes, there are two other 
virtual modes that involve the field quanta with negative 
frequencies, i.e., u>' v = — (u>i+u)2) and u v = — (wi — W2) < 
0. They are initiated by the interaction terms a\a2 a k + 
h.c. and a\a 2 ak + h.c. of Eq.(Bl), respectively. The u' v 
mode is a fast rotating mode so that its contribution to 
the reduced dynamics averages to zero in the relaxation 
time scale, tr <~ I/7, where the A 2 1- approximation is 
valid. This is the rotating-wave approximation [7, 40] 
usually implemented in the Markovian limit. 

The virtual mode u v does not contribute to the re- 
duced dynamics on the level of the A 2 i-approximation 
with the collision operator ip 2 defined in Eq.(Cl). A 
similar example is provided by Ref. 23 for the reduced 
dynamics of a single oscillator coupled to a field, in which 
the interactions contain a virtual transition mode. The 
effect of a possible extension of the collision operator ^ 
(CI) to include the contribution of the virtual transitions 
is presented in the next appendix. 

It is also interesting to note that if instead of assuming 
lj\ > LJ2 in our discussion so far, the opposite situation 
L02 > to 1 would interchange the role played by the oj v 
and the wo modes, i.e., now ujq becomes a virtual mode, 
whereas uo v becomes a resonant mode. However, the uj' 
and uj' v modes are not affected by this change. 

In summary, in comparison to the full anharmonic in- 
teractions, the higher frequency resonant mode u' is not 
included in the interaction Hamiltonian of the system we 
consider here. The contribution of this mode is negligible 
compared to the loq mode at low temperature, or when 
a frequency cut-off is imposed to the form factor below 
cu' . Moreover, out of the two virtual modes we have dis- 
cussed, the effect of the uj' v mode is effectively dropped 
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FIG. 3. p("i."2)^ 2 p("i -2,^2+2) transition . It is one f the 
diagrams giving the contribution to the collision operator %j)2 
when virtual mode w„ is considered. The operators on the 
right vertex cause the virtual transition. 



under the rotating-wave approximation, whereas the ui v 
mode does not affect the dissipative part of the reduced 
dynamics using the usual definition of the collision oper- 
ator. 



Appendix C: Extension of ip2 

Using the usual definition of the collision operator up 
to the second order in the coupling constant [23], 



r 2 = p v ^p v 



P P C P" - \ 2 P"£i 



1 



£0 — w • v 



(CI) 

fc v p° , 

(C2) 



the virtual mode u> v does not contribute to the collision 
operator on the level of A 2 i-approximation. However, if 
we incorporate the transitions between different P v sub- 
spaces into the definition of the collision operator (CI), 
the virtual mode u v may contribute to the dissipative 
part of the reduced dynamics as follows [41]. 

We define the collision operator by i\j' 2 = PipzP, where 
P = ^2jpP v ■ Fig. 3 shows an example of a possible tran- 
sition of -02 that involves different P v subspaces. As indi- 
cated in the figure, the interaction vertex on the right is 
due to the virtual transition a[a2a k , whereas the vertex 
on the left is due to a real transition mode in the origi- 
nal model (2). A standard calculation shows that Fig. 3 
contributes a term (a}o2) 2 / to the dissipative part of the 
reduced dynamics, in addition to the contributions by 
the real transitions. When taking into account all these 
additional transitions, the terms L + fL + , L 2 + f , jP\ and 
their hermitian conjugates, are added to the dissipative 
operator, resulting in a not completely positive reduced 
dynamics [42]. 

We note that this extension to the definition of the 
collision operator does not affect the reduced dynamics 
of the original model (2), since all the possible transitions 
in this model involve the same P v subspace. We already 
learned that in this situation they lead to a completely 
positive reduced dynamics [26, 27] with the Kossakowski- 
Lindblad form (9). 
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Appendix D: Rotational symmetry under Lq 



In this appendix, we show that K is invariant under a 
rotation along the Lq axis. Indeed, we find that 



L' ± = e i6Lo L ± e- i$L ° = e ±t8 L± , 



(Dl) 



the /( JV > jV ) basis, we have 



.v 



r n{t 



(N,N) r(N,N) 

" 2n;N-2m ' 



= STA N ' N W- N - 

/ j °w;m JN — 



^ ) = ^( Sln2 </©( 



ml 



(E2a) 
(E2b) 



by using the commutation relations [L + ,LJ\ = 2Lq and 
[Lq,L±] — ±L±. Since L+ and _L_ appear pairwise in 
K, cf. (6)-(9), the phase in Eq.(Dl) cancels out. Added 
by the fact that Lq commutes with N, we conclude that 
K is invariant under the rotation. 



Appendix E: SU(2) generalized coherent states 
reside in /' ' ' subspace 



The SU(2) generalized coherent states [30-32] 



N 

\r) N = (l + \r\ 2 )- f Yl 
m=0 




T ni \m,N-ni). 



(El) 



reside in the corresponding f( N < N > subspace of the sys- 
tem, where r = tan(0/2) cxp(— i<p), in which 0, </> are the 
parameters that label the coherent states. In terms of 



Appendix F: Kd in the interaction picture 

In this appendix, we calculate the effect of Kd on the 
fi basis. Since Lq and N commute, we have 

ift'gt _ i(u:' -Su:' )Nt/2-iSu:' N 2 t/4 iSu)' g Llt i(uo-Sui )L t 

(Fl) 

Using Eq.(Dl), we find that 

Therefore, K d = e^o^/Qe^'K^o*, since L + and L_ 
appear pairwise in Kd, and N commutes with Lq,L±. 
Furthermore, using 

[Ll, L ± ] = ±(2L T l)L± = ±L ± (2L ± 1) , (F3) 

we can show that 

L± = e lSu '' aL » t L±e~ iSuj '° L » t = 7 J±e ±* <5w o( 2L o±i)t rp4\ 

= p ±iS^' (2L 0T l)t L± ^ ( F5 ) 

where L* ± = L T . Using the relation (17c), we finally 
obtain Eq.(39). Therefore, Kd differs from Kd by a phase 
factor in two of the terms. For basis vectors that lie in 
the probability subspace, v = r — f = 0, and we get 
K d = K d . 
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